New code for equilibriums and quasiequilibrium initial data of compact objects. II. 
Convergence tests and comparisons of binary black hole initial data 
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COCAL is a code for computing equilibriums or quasiequilibrium initial data of single or binary 
compact objects based on finite difference methods. We present the results of supplementary conver- 
gence tests of COCAL code using time symmetric binary black hole data (Brill-Lindquist solution). 
Then, we compare the initial data of binary black holes on the conformally flat spatial slice obtained 
from COCAL and KADATH, where KADATH is a library for solving a wide class of problems in theoretical 
physics including relativistic compact objects with spectral methods. Data calculated from the two 
codes converge nicely towards each other, for close as well as largely separated circular orbits of 
binary black holes. Finally, as an example, a sequence of equal mass binary black hole initial data 
with corotating spins is calculated and compared with data in the literature. 



I. INTRODUCTION 

Various methods have been developed in the past cou- 
ple of decades for computing numerical solutions of com- 
pact objects in equilibrium or quasiequilibrium. Those 
include methods for computing relativistic rotating stars 
in equilibrium or binary black hole initial data (see e.g. 
[J, [J). For example, numerical solutions of binary black 
holes (BBH) in quasicircular orbits have been widely used 
for initial data of merger simulations and sequences 

of such data with fixed irreducible mass of each black hole 
(BH) have been also applied to approximate the inspiral 
evolution of BBH [|. 

The authors have been developing, independently, nu- 
merical codes for computing such compact objects, CO- 
CAL [i| (Paper I hereafter) and KADATH (Paper II). 
CoCAL is a code for computing various kinds of astro- 
physical compact objects - isolated or binary systems 
of neutron stars and black holes which may be associ- 
ated with strong magnetic fields. KADATH is a library for 
solving a wide class of problems in theoretical physics 
including those of general relativity, and is capable of 
computing such compact objects. 

In the first part of this paper, we present the results of 
supplementary convergence tests of COCAL to those pre- 
sented in Paper I. With straightforward changes in the 
radial coordinate grid spacings and in the finite differ- 
ence formula for the integration over the zenith angle, 
errors in the gravitational fields especially near the com- 
pact objects decrease substantially, which is necessary to 
improve the accuracy of widely separated BBH solutions. 
In the second part of the paper, we carefully compare the 
spatially conformally flat BBH initial data in circular or- 



bit calculated from COCAL and KADATH code. Compari- 
son of the solutions is the most effective test to confirm 
the reliability of the codes in which a system of com- 
plicated equations is solved. Such comparison had been 
done for the codes for relativistic rotating stars in 
To our knowledge, this is the first attempt to compare 
the BBH data calculated from totally different methods 
as the spectral method and the finite difference method. 
Finally, we present a sequence of spatially conformally 
flat BBH initial data in circular orbits for the case with 
equal mass and corotating spin, and compare the result 
with those presented in [5|. Throughout the paper we 
use geometric units with G = c = 1. 



II. CONVERGENCE TESTS FOR COCAL CODE 

In this section, we present convergence tests of COCAL 
code supplementary to those presented in Paper I. The 
setup for the test problem is the same as in Paper I: 
the Brill-Lindquist solution for the time symmetric BBH 
data is generated numerically, and it is compared with 
the analytic value. We briefly review the setup of the 
problem and discuss the modification of the finite differ- 
ence scheme from the previous paper. 



A. A test problem for binary black holes 

We assume the spacetime Ai is foliated by a fam- 
ily of spacelike hypersurfaces (Et)tgR, Ai = R x Y, 
parametrized by t g R. We assume the line element 
in the neighborhood of Sf to be 
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where fij is the flat spatial metric, so that the data on 
Yit becomes time symmetric - the extrinsic curvature Kab 
on Ei vanishes. 



2 



Decomposing Einstein's equation Gap = with re- 
spect to the fohation using hypersurface normal to 
St, and the projection tensor 7°'' g"^ + n^n^ to it, 
we write the Hamiltonian constraint Gapn"n^ = 0, and 
a combination of the spatial trace of Einstein's equation 
and the constraint Gapi'^"^ + ^n"n^) = 0, as 



= and V2(aV) = , 



(2) 



These equations have solutions, which correspond to the 
Schwarzschild metric in isotropic coordinates for a single 
BH. For a two BH case, a BBH solution is given by Brill 
and Lindquist Q: 

Ml M2 , , Ml M2 , , 

^- = 1 + -^ + -^ and ct/. = 1 - - 3 
2ri 2r2 2ri 2r2 

where subscripts 1 and 2 corresponds to those of the first 
and second BH; ri and r2 are distances from the first 
and second BH, respectively, and Mi and M2 are mass 
parameters. Instead of solving two Laplace equations 
Eq. (U), we write an equation for a with a source on the 
whole domain of Et: 



= and 



(4) 



In an actual computation, spherical regions near the cen- 
ter of BH are excised to avoid singularities. Therefore, 
boundary conditions for these elliptic equations are im- 
posed at the radius r = of the excised sphere Sa, and 
at the radius r = of the boundary of computational 
domain Sb- We also set the mass parameters Mi and M2 
as 0.8 X Tq of each BH to avoid the lapse to be negative at 
Sa- In the following tests, we impose Dirichlet boundary 
conditions on Sa and Sb whose values are taken from the 
analytic solution 



B. Coordinates, grid setup, and finite difference 
scheme of COCAL code 



As explained in Paper I, three spherical coordinate 
patches are introduced for solving binary compact ob- 
jects with COCAL . Two of them arc the compact ob- 
jects coordinate patches (COCP-I and II) and one is 
the asymptotic region coordinate patch (ARCP). In each 
spherical patch, coordinates cover the region (r, 9, (f>) G 
['^aj'^b] X [0,7r] X [0,27r]. The two COCPs are centered at 
the center of compact objects and extend up to about 
rb ~ 0(10'^ M), while ARCP is centered at the center of 
mass of the binary, and extends from ^ O(10il/) to 
Tb ^ O(IO^M), where M is the total mass of the sys- 
tem. Definitions of the parameters for the grid setups 
arc listed in Table HI 

In solving a system of elliptic equations such as Eq. dH , 
we rewrite them in integral form using Green's function 
that satisfies given boundary conditions, and apply a fi- 
nite difference scheme to discrctize those integral equa- 
tions on the spherical coordinates of each domain. We use 



Va : Radial coordinate where the radial grids start. 
r(, : Radial coordinate where the radial grids end. 
Tc : Radial coordinate between Va and rt where 

the radial grid spacing changes. 
re : Radius of the excised sphere. 
Nr : Number of intervals Ar^ in r G [ra , r-;,]. 
NI : Number of intervals Art in r G [ra, 1]. 

: Number of intervals Avi in r G [va , r^] . 
Ne : Number of intervals A9j in 9 e [0, tt] . 

: Number of intervals A(j}k in G [0, 2tt]. 
d : Coordinate distance between the center of Sa (r — 0) 

and the center of mass. 
da : Coordinate distance between the center of Sa {r — 0) 

and the center of Se- 
L : Order of included multipoles. 

TABLE I. Summary of grid parameters for COCP. 



the midpoint rule for numerical quadrature formula, and 
hence compute the source terms at the midpoints of the 
grids. Coordinate grids (ri,6j,4>k) with i = 0, • • • ,iVr, 
j = 0, ■ ■ ■ , Ng, and fc = 0, • • • , A^^, are freely specifiable 
except for the endpoint of each coordinate grid that cor- 
responds to the boundary of the computational region, 
(ro, 6*0,00) = (fa, 0,0) and {rN,,0Ne,4>Nj = (^'fc, ti", Stt). 
The grid setup for COCP and ARCP is the same as Pa- 
per I except for the radial grid of COCP which will be 
explained later. For angular coordinate grids {9j, 0^), we 
choose equally spaced grids. 

In Paper I, we have used for the finite difference 
formulas, (1) 2nd order midpoint rule for the quadra- 
ture formula, (2) 2nd order finite difference formula for 
the 9 and (j) derivatives evaluated at the mid points 
(r,+ i,0^-+i,0fc+i) = ((r, +ri+i)/2,(0i + 0,+i)/2,(0, + 
0i+i)/2), (3) 3rd order finite difference formula for the 
r derivative evaluated at the mid points, and (4) 4th or- 
der finite difference formula for the derivatives evaluated 
at the grid points, {ri,6j,4>k)- In the present compu- 
tations, we use the same finite difference formulas men- 
tioned above except for the numerical quadrature formula 
in 9 integrations. 

Differences from the previous Paper I are the spacings 
of radial grids Ar^ := — r^-i in COCP, and the quadra- 
ture formula used for the integration in zenith angle 9. 



1. Radial grid spacings for COCP 

When we compute a sequence of BBH data from larger 
to smaller separations in the COCAL code, we change the 
BH excision radius from smaller to larger values and 
fix the separation ds (instead of fixing and varying 
dg). In this way, the number of grid points are kept to 
be the same, and the structures of coordinate grids are 
almost the same for all solutions of the sequence. As 
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FIG. 1. The radial grid points ri of COCP are plotted against 
the grid number i = 0, • • • , Nr for the case with H3 grid in 
Table El 
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Values of the parameters for the coordinate grids of 
COCAL used in computing the results presented in this 
paper are listed in Table HIl In Fig. [1] an example of the 
radial grid points is plotted for the case with 113 grid 
setup in Table HT] Because of the construction, the grid 
structure in the region larger than r > 1 is the same for 
all solutions with different BH radius once a grid setup 
(resolution) as in Table HIl is selected. 



a result the discretization error behaves systematically 
from one solution to the other, and hence the quantities 
such as mass or angular momentum vary smoothly along 
a sequence of solutions. 

It is important to notice that, by changing the BH 
radius, we change the mass of the solution and hence 
change the length scale of the system. Therefore, to 
maintain the accuracy of the gravitational fields near the 
BH, the intervals near the hole should be proportional to 
the mass of the BH, or in our case, the BH excision ra- 
dius Ta- Therefore, we modify the construction of the grid 
spacing in the radial direction r of COCP for computing 
a sequence from smaller to larger BH as follows. Without 
loss of generality, we set the radius of BH excision sphere 
Sa as Ta < 1. We divide the radial coordinate to four 
regions, I: r G [va, 1], II: r G [l,rc], HI: r G [rc,3rc], and 
IV: r G [3rc,r(,]. We set the first interval by 

An ^ ^, (5) 

where is the number of intervals in the region I: [r^ , 1] , 
and A is a constant factor which is chosen to be A = 0.75. 
For each region, Ari fi — ^i-i, are defined by 

An+i = hiAn, for i = 1, ■ • • , iV^^ - 1 (6) 

Ar, = Ar, for i^N^,--- ,iV;" (7) 

An+i = h3An, for i ^ N';\- ■ ■ + - 1 (8) 

Ar,+i = hiAn, for i = iV^" + iV^ , • • • , iV, - 1 (9) 

which correspond to regions I, II, III, and IV, respec- 
tively. The ratios hi{> 1) (i = 1,3,4) are respectively 



2. 4th order midpoint rule for the quadrature formula of 9 
integration 

As discussed in Paper I, our Poisson solver is a sys- 
tem of integral equations, and it is numerically integrated 
with a quadrature formula of midpoint rule. Therefore, 
the sources of the integrals are always evaluated at the 
midpoints of (r.^, 0j, (j)k) grids. As summarized above, the 
2nd order midpoint rule was used for a quadrature for- 
mula in Paper I. With the above mentioned choice for 
finite difference formulas, the 2nd order convergence of 
the error has been achieved. As shown in the top panel 
of Fig. |2l (as well as figures 3, 4, and 6 in Paper I), how- 
ever, the fractional error of the potentials normally in- 
creases near the excision surfaces of the BH, Sa {r = ra), 
although it converges in 2nd order. 

One might expect that the increase of the error near Sa 
is due only to a lack of resolution in radial grid points. It 
turns out, however, that the finite difference errors in the 
potentials near the boundaries of computational domains 
arc dominated by the discretization error in the 6 coor- 
dinate. In particular, the 6 integration of the source in- 
volving the Legendre function turns out to be the source 
of error. Therefore, we replace the quadrature formula of 
9 integration to 4th order accurate midpoint rule whose 
weights are 




where the grid number j is a multiple of 4, and AO = 
TT/Ne. 
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Type 


Patch 


Ta 




Tc 


Te 


Nl 






Ng 




L 


HI 


COCP-1 


var. 


10^ 


1.25 


1.125 


32 


40 


96 


24 


24 


12 




COCP-2 


var. 


10^ 


1.25 


1.125 


32 


40 


96 


24 


24 


12 




ARCP 


5.0 


10^ 


6.25 


— 


8 


10 


96 


24 


24 


12 


H2 


COCP-1 


var. 


10^ 


1.25 


1.125 


64 


80 


192 


48 


48 


12 




COCP-2 


var. 


10^ 


1.25 


1.125 


64 


80 


192 


48 


48 


12 




ARCP 


5.0 


10^ 


6.25 


— 


16 


20 


192 


48 


48 


12 


H3 


COCP-1 


var. 


10^ 


1.25 


1.125 


128 


160 


384 


96 


96 


12 




COCP-2 


var. 


10^ 


1.25 


1.125 


128 


160 


384 


96 


96 


12 




ARCP 


5.0 


10^ 


6.25 




32 


40 


384 


96 


96 


12 


H4 


COCP-1 


var. 


10^ 


1.25 


1.125 


256 


320 


768 


192 


192 


12 




COCP-2 


var. 


10^ 


1.25 


1.125 


256 


320 


768 


192 


192 


12 




ARCP 


5.0 


10" 


6.25 




64 


80 


768 


192 


192 


12 



TABLE II. Grid parameters of COCAL used for computation 
of BBH data. The separation of two BHs is fixed as ds — 2.5. 
For the excision radius for COCP-1 and 2, "var." stands 
for a variable parameter assigned to each solution. In the test 
problems in Sec lII CI they are chosen to be ra = 0.2 and 0.4 
for close BBH, and = 0.05 and 0.1 for separated BBH for 
COCP-1 and 2, respectively. 



C. Convergence tests 

We perform convergence tests to examine that the 
above two modifications improve the accuracy of the CO- 
CAL code. In Table [ITl grid setups for the computations 
are listed. The grids H1-H4 correspond to different levels 
of resolutions. At each level, the resolution is double the 
previous onc0. 

In Figs. [2H11 the fractional errors in the lapse that 
are averaged over the angular coordinate grids {9i^(j)j) 
at fixed radial coordinate r are plotted against r of each 
coordinate patch, 



5a 
a 



1 



E 



^ Ci; exact 



Q^cxact 



(14) 



where writing a grid point {ri,9j,(j)k) by p, we define a 
set Qi by Qi := {p \ p S™ and = const}, where 

5™ is an interior domain of Se, and ^{Qi) is the number 
of points included in Qi. 

In the top panel of Figl^J the same finite difference 
scheme as presented in Paper I is used for computing 
closer BBH solution with the BH excision radius = 0.2 
and 0.4, and with the separation ds = 2.5. The midpoint 
rule in the 6 integration is 2nd order accurate in this 
panel. In the bottom panel, the 4th order midpoint rule 
Eq. ([T3| for the 9 integration is used for the same model 



We have also tested different combinations of grid numbers 
(Nr,Ng, N^) for the first ievel of resolution and performed con- 
vergence tests. We found the combination of type H was better 
than others. For example, the accuracy was not improved by 
increasing the grid points in the (p direction. 
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FIG. 2. Angular averaged fractional errors in the lapse 
{|5a/a|) are plotted along the radial coordinate r for the 
nonequal mass BBH data computed on three multiple 
patches. Top panel: data computed with the same coordi- 
nate grid spacings and finite difference schemes as in Paper I. 
Bottom panel: data computed with the 4th order 6 integra- 
tion discussed in Sec lII B 2l The grid parameters and number 
of grid points are varied as H1-H4 in Table ITIl The BH exci- 
sion radii are chosen to be ra ~ 0.2 and 0.4 for COCP-1 and 
2, respectively, and the separation is to be da = 2.5. 



and the same grid spacings. Clearly, the fractional error 
substantially decreases by this change for the H2 to H4 
levels. We notice that the error near the BH converges 
in 4th order, that is, the error decreases about 1/16 at 
each level of resolution. The errors near the BH as well 
as in the asymptotic region are dominated by those from 
the 9 integrations of the surface integral terms. 

In FiglSl we calculate more separated BBH solutions, 
decreasing the radius of BH excision surface to 1/4 
(therefore effectively separating BBH 4 times apart) as 
Va = 0.05 and 0.1 with the same separation dg = 2.5, 
and using the 4th order integration in 9 as in the bottom 
panel of Figl^] Although the errors near the BH excision 
surfaces are of the same order of magnitude as those of 
the corresponding resolutions plotted in Fig. [21 bottom 
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FIG. 3. Same as Fig. [g] but for the BH radius ra = 0.05 and 
0.1 for COCP-1 and 2, respectively. The resuhs are calculated 
using the 4th order midpoint rule for integration as in the 
bottom panel of Fig[2l otherwise the same finite differencing 
scheme (in particular, the same radial grid spacing Ar^) as in 
Paper I. 



panel, the errors once increase as the radial coordinate 
r increases. It turns out that the BBH initial data dis- 
cussed later can not be calculated accurately with this 
grid setup for largely separated orbits. 

In FiglH convergence tests for the close {va = 0.2 and 
0.4) and the separated (r^ = 0.05 and 0.1) BBH are 
calculated with scaled radial spacing near the BH intro- 
duced in Scc lIIB ll as well as the 4th order midpoint 
rule in 9 integration. The size of the errors around the 
BH excision radius for the largely separated BBH data 
(bottom panel) is now comparable to those for the close 
BBH data (top panel) for each level of resolution. This 
improvement turns out to be important for accurately 
computing the separated BBH data in COCAL code. 
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FIG. 4. Same as Fig. (2] The results are calculated using 4th 
order integration in 6 coordinate Eq. (|13[1. and the scaled ra- 
dial spacing discussed in Sec lII B T] Top panel: data with the 
BH radius ra = 0.2 and 0.4 for COCP-1 and 2, respectively. 
Bottom panel: data with the BH radius = 0.05 and 0.1 for 
COCP-1 and 2, respectively. 



III. COMPARISON OF BBH INITIAL DATA 

A. The KADATH library 

In this section, wc compare the circular solutions of 
BBH initial data on a conformally flat spacelike hyper- 
surface calculated from COCAL and the KADATH library. 
KADATH 01 is a library designed to solve a wide class of 
problems in theoretical physics including those of general 
relativity such as the above compact objects. It is based 
on spectral methods (sec for instance [10| and references 
therein) where the various fields are approximated by fi- 
nite sums of known functions typically trigonometrical 
functions and orthogonal polynomials. One of the main 
advantages of spectral methods is their fast convergence 
to the true solution (typically exponentially), when one 
increases the order of the expansion. For instance, in this 
paper, a relative accuracy of about 10~^ is achieved with 



15 coefficients in each dimension. 

Spectral methods enable one to translate a set of par- 
tial differential equations into an algebraic system on the 
coefficients of the expansions. This system is then solved 
by a standard Newton-Raphson iteration. The compu- 
tation of the Jacobian as well as its inversion are paral- 
lelized. 

The code used in this paper is essentially the same as 
the one used in Sec. 7.3 of In order to check the 
overall accuracy of the computations, one monitors the 
convergence of some global quantities (like the orbital fre- 
quency), as a function of N, the number of points in each 
dimensions. Let us mention that, in the case of a large 
separation, the code was slightly modified to maintain 
accuracy, probably due to a stretch of the bispherical co- 
ordinates when the distance between the holes gets much 
bigger than the size of the holes themselves. In particu- 
lar, a spherical shell was added between the bispherical 
coordinates and the outer compactified domain and the 
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determination of the orbital velocity had to be changed 
(see Sec. HITC]) . 



B. Conformally flat BBH initial data 

The circular solution of BBH initial data is calculated 
by solving the Hamiltonian and momentum constraints, 
and the spatial trace of the Einstein's equation on a con- 
formally flat spacelike hypersurface Et. The spacetime 
metric on St is written in 3+1 form as 



Code 



Res. QMir 



MADM/A/irr J/Mj^rr Mjrr/rq 



ds = g^^dx'^dx" 

= -a^dt^ + Jijidx'' 



- I3'dt){dx' + P^dt), (15) 



where the spatial three metric 7^ on the slice St is as- 
sumed to be 7y = ip'^fij. Here, field variables ^^tU^ and 
/3' are the conformal factor, lapse, and shift vector, re- 
spectively, and fij is a fiat three dimensional metric. We 
also assume maximal slicing to Et, so that the trace of the 
extrinsic curvature Kij := — ^(£t7y — ^plij) vanishes. 
Writing its tracefree part Aij, the conformally rescaled 
quantity Aij becomes 



Aij — 



1 

2^ 



(16) 



where the derivative di is associated with the flat metric 
fij, and conformally rescaled quantities with tilde are 
defined by Ai^ — A.i^ and = ^% whose indexes are 
lowered (raised) by fij (/'^). The system to be solved, 
which are Hamiltonian and momentum constraints and 
the spatial trace of the Einstein's equation, becomes 



At/; 



Aij A^-' , 



= -2aA,^d, In ^ - ^5,5, 
a 3 

A(a^) = laij'A,jA'\ 

o 



(17) 
(18) 
(19) 



where A := 9^9* is a flat Laplacian [TTI. [l^ . 

For the boundary conditions at the BH excision bound- 
ary Sa, we choose approximate irrotational apparent 
horizon boundary conditions. 



9-0 
dr 



2r 



^3 



no, 



(20) 

(21) 
(22) 



where no is an arbitrary positive constant for which we 
choose no = 0.1, is the unit normal to the sphere Sa, 
and fl represents a parameter for orbital angular veloc- 
ity. The vector yl^^ := (0, d, 0) is the translational vector 
with respect to the center of mass. With these condi- 
tions, the sphere Sa becomes an apparent horizon (AH) 
in quasiequilibrium js], 0, [l^ . 



KADATH 
KADATH 
KADATH 

COCAL 
COCAL 
error [%1 



11 0.127171 

13 0.127299 

15 0.127346 

H2 0.127243 

H3 0.127383 

- 0.03 



ds/ra = 12 
0.982866 
0.983041 
0.983072 
0.982848 
0.983125 

0.005 



0.757608 3.93366 

0.758390 3.93405 

0.758698 3.93417 

0.754592 3.93191 

0.758389 3.93392 

0.04 0.006 



KADATH 
KADATH 
KADATH 

COCAL 
COCAL 
error \% 



11 0.0433696 

13 0.0350000 

15 0.0349704 

H2 0.0344185 

H3 0.0348624 

- 0.3 



ds/Va = 30 

0.933609 
0.988279 
0.988129 
0.987369 
0.987982 
0.01 



1.129241 3.56317 

0.906819 3.55178 

0.907249 3.55166 

0.889918 3.54574 

0.903226 3.55021 

0.4 0.04 



TABLE 111. Comparison of BBH initial data from COCAL and 
KADATH for the cases ds/va = 12 and ds/ra = 30. Fractional 
differences in % between the highest resolution cases of KADATH 
and COCAL are shown in the lines indicated by "error". A 
column "Res." stands for the resolution. 



At the asymptotics, the boundary conditions are 





.00 - 1-0 




= 0.0 

-foo 







(23) 
(24) 
(25) 



When using KADATH , the whole spacelike slice St is com- 
pactified, and hence the above conditions are imposed at 
the spatial infinity, while in COCAL , the computational 
domain is truncated at the radius ^ O(IO^M) and it 
is at rjy that the above conditions are imposed. 



Comparison of the circular solutions of BBH 
initial data for KADATH and COCAL 



Following [IJ], we obtain the angular velocity f2 of 
a circular orbit of BBH initial data from an assump- 
tion that an equality of ADM mass and Komar mass, 
A^ADM = A/k, is satisfied for the circular orbit, where 
AfADM and Mk are defined by 



A/a 



1 j ^jacfM _ fabfcd-^Q^^^^dSa 

1671- Jsoc 

- / D^'HSa, (26) 
Mk = -^ f V^t^dSc.0 

D^'adSa- (27) 



1 

An 



In COCAL , the surface integrals are calculated at a cer- 
tain finite radius, typically around r (10^ Af). 
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FIG. 5. The panels show the scaled values of MADM/Mirr as 
a function of Mirr/fa and Jadm/Mj^j as a function of QAfiir- 
The circles denote the results from KADATH and the triangles 
the ones from COCAL . The runs using KADATH are labeled 
by the number of points in each dimension. The left top and 
bottom panels correspond to a separation of ds/va ~ 12 and 
the right ones to ds/va ~ 30 (in this case the lower resolution 
results from KADATH are outside the range of the plot). 



In KADATH , these integrals arc usually evaluated at 
spatial infinity as its definition, / :~ lim / with 

Sr the sphere of radius r. However, in the case of a 
large separation, it turned out that this was not giving 
good results. The precision of the results, measured by 
convergence of the value of fi, is much better when one 
demands that, at spatial infinity. 



1 



0(r- 



(28) 



The difference between the Komar and ADM mass is then 
of the order of 10~^ for the higher resolution, thus giving 
a measure of the overall error of the code. 

For a converged circular solution, we also calculate 
ADM angular momentum: 



J 



1 

8^ 



(29) 



FIG. 6. Fractional errors in metric potentials between those 
calculated by KADATH and COCAL . Errors in the conformal 
factor ^, lapse q, and j/-component of the shift Py are plotted 
along the s-axis that intersects with the centers of two BHs. 
The model is the same as that in Table HIIl with the separation 
ds/ra = 12. 



The above quantities are normalized by the irreducible 
mass of AH, Mirr, which is defined from the surface area 
of AH, namely the area integral over the BH excision 
surface Sa, 



(30) 



We write Mhri := -\/^AHi/167r and Afiir2 '■= 
^/]4AH2/167r for each BH, and write A/in- = + Miy^2 
for a total mass. 

In Table Hill global quantities normalized by Mirr are 
presented for the irrotational BBH data computed from 
KADATH and COCAL at the separations ds/ra = 12 and 
ds/ra = 30. Three different resolutions for KADATH are 
given, mainly 11, 13 and 15 points in each dimension, 
and two (lower and higher) resolutions, H2 and H3, are 
used for the computations of COCAL . The relative dif- 
ferences between the highest resolution results from both 
codes are also indicated. The convergence of these quan- 
tities is plotted in Fig.[5]for both separations, dg/ra = 12 
(left panels) and ds/ra = 30 (right panels). All plots in- 
dicate a nice convergence of the global quantities, when 
the resolutions are increased. 

In FiglH plotted are fractional errors of metric poten- 
tials of the same BBH data between those calculated 
from KADATH and COCAL . The errors are defined by 
100 X |g(C0CAL) - g(KADATH)|/|g(KADATH)| [%], and those 
of the conformal factor ip, lapse a, and ?/-component of 
the shift Py are plotted along the a;-axis which intersects 
with the centers of two BHs for the case with the sep- 
aration ds/ra ~ 12. Resolutions are H3 for COCAL and 
15 points for KADATH . As seen from the figure, the met- 
ric potentials from the two codes agrees well. Note that 
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FIG. 7. Sequence of solutions for conformally flat BBH initial 
data with corotating spins are plotted. The sequence com- 
puted from COCAL code is compared with those presented in 
Caudill et.al. [a| (CCGP) and with the third post-Newtonian 
(3PN) results 



the relatively large error in (3y near x/ra « 6 is due to 
(3y crossing zero, and hence its fractional error diverges 
there. 



D. Solution sequence for corotating BBH data 

Finally, we present a sequence of solutions for the con- 
formally flat BBH initial data computed from COCAL . In 
this computation, the boundary condition of shift /?* (|2T|) 
is replaced by that for the BBH with corotating spins. 



/3l 

' r — r^ 



"0 » 



(31) 



where the vector (f>l^ is a generator of rotation around 
the center of mass whose components in Cartesian co- 
ordinates is written (pl^ = (— j/cm, a;cm, 0). The orbital 
angular velocity parameter n is evaluated with the same 
method as mentioned in the previous section. 



In Fig. [3 a sequence of the corotating BBH data com- 
puted from COCAL with resolution H3 is compared with 
those of the paper by Caudill et.al. [1], and of the third 
post-Newtonian (3PN) circular orbits [l^. In the top 
panel, the binding energy normalized by the irreducible 
mass Eb/Mi^y := (A/adm — Mirr)/Mirr is plotted against 
the normalized angular velocity fiAfii-r, and similarly in 
the bottom panel, the angular momentum J/Mf^^. The 
curves calculated from COCAL agree well with the other 
two curves. In the smaller flMiir (large separation), the 
curves from the COCAL are slightly smaller than those 
of the other works. The size of this error is comparable 
to that listed in Table HhI The error in Ei,/Mi^^ for the 
model with dg/ra = 30 is around —1.2% whose separa- 
tion corresponds to the one in Figl7]with VlMi„ ~ 0.026 
{ds/va — 31.25). The data used in Fig. [7] are tabulated 
in Table. |lVl 

For computing a solution with the COCAL code, each 
iteration takes about 3 minutes for H3 resolution with a 
single CPU (1 core) of Xeon X5690 3.46GHz, and each 
run uses about 6GB of RAM. A convergence to a circular 
solution is achieved after 500-700 iterations, during which 
an iterative search for the circular f2 is made 5-7 times. 



IV. DISCUSSION 

We have presented additional convergence tests of CO- 
CAL code focusing on the BBH data. Especially the con- 
formally flat initial data of BBH in circular orbit calcu- 
lated from COCAL code are compared with those from 
KADATH code. We demonstrated that the results from 
both codes converge toward each other for large and small 
separations (see also [IB])- 

As fully discussed in Q, a corotating sequence 
presented in Sec lIIIDl is not considered model for 
inspiraling BBH, because BBH tides do not effectively 
work to spin up the BH to synchronize the BH spin 
with orbital motion within the time of inspirals. In [5|, 
the authors describe a more realistic sequence of BBH 
inspiral where the spin angular momentum of the AH is 
conserved. We will present elsewhere the performance 
of the COCAL code for computing those sequences to 
model BBH inspirals, which is a necessary step to 
compute more complex binary systems including black 
hole-neutron star binaries. 
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Ta ds/Va 

0.03 83.333 

0.04 62.500 

0.05 50.000 

0.06 41.667 

0.07 35.714 

0.08 31.250 

0.09 27.778 

0.10 25.000 

0.11 22.727 

0.12 20.833 

0.13 19.231 

0.14 17.857 

0.15 16.667 

0.16 15.625 

0.17 14.706 

0.18 13.889 

0.19 13.158 

0.20 12.500 



24.327 
18.126 
14.403 
11.921 
10.147 
8.8151 
7.7785 
6.9483 
6.2682 
5.7007 
5.2197 
4.8067 
4.4481 
4.1337 
3.8557 
3.6081 
3.3860 
3.1856 



0.0077690 
0.011925 
0.016575 
0.021647 
0.027086 
0.032848 
0.038898 
0.045207 
0.051746 
0.058495 
0.065430 
0.072533 
0.079787 
0.087174 
0.094679 
0.10229 
0.10998 
0.11775 



Madm/M r 
0.99501 
0.99360 
0.99226 
0.99100 
0.98984 
0.98877 
0.98780 
0.98693 
0.98615 
0.98548 
0.98491 
0.98445 
0.98408 
0.98381 
0.98364 
0.98356 
0.98358 
0.98368 



Eb/Mirr 

-0.0049873 
-0.0064032 
-0.0077416 
-0.0089952 
-0.010159 
-0.011228 
-0.012201 
-0.013074 
-0.013847 
-0.014518 
-0.015087 
-0.015555 
-0.015922 
-0.016191 
-0.016362 
-0.016439 
-0.016424 
-0.016320 



1.3211 
1.1824 
1.0910 
1.0268 
0.97988 
0.94470 
0.91792 
0.89736 
0.88156 
0.86946 
0.86032 
0.85355 
0.84873 
0.84550 
0.84360 
0.84279 
0.84290 
0.84377 



0.10277 
0.13793 
0.17357 
0.20972 
0.24638 
0.28360 
0.32140 
0.35980 
0.39884 
0.43854 
0.47896 
0.52011 
0.56204 
0.60479 
0.64839 
0.69289 
0.73834 
0.78477 



TABLE IV. Solution sequence of equal mass BBH with corotating spins. In the computation, parameters H3 listed in Table 
in] are used with a fixed separation ds = 2.5, and with varying the excision radius Va of BH. 
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